library("lfe")
library(texreg)
library(readstata13)
library(plm)
library(Hmisc)
library(haven)
library(dplyr)
library(plyr)



setwd("C:\\Users\\Public\\Documents\\Results\\Establishment_Composition")

simplelag <- function (x,by=NULL,mylag=1,outside=NA) {
  myend <- length(x)- mylag
  if(!is.null(by))  {
    lby <- c(replicate(mylag,""),as.character(by[1:myend]))
    y0 <- c(replicate(mylag,outside),x[1:myend])
    y <- ifelse(as.character(by)==lby,y0,outside)
  }
  else {
    y <- c(replicate(mylag,outside),x[1:myend])
  }
}

forward <- function (x,by=NULL,myforward=1,outside=NA) {
  mystart <- myforward+1
  if(!is.null(by))  {
    fby <- c(as.character(by[mystart:length(x)],replicate(myforward,"")))
    y0 <- c(x[mystart:length(x)],replicate(myforward,outside))
    y <- ifelse(as.character(by)==fby,y0,outside)
  }
  else {
    y <- c(x[mystart:length(x)],replicate(myforward,outside))
  }
}


retain <- function(x,event,outside=NA)
{
  indices <- c(1,which(event==TRUE),length(x)+1)
  values <- c(outside,x[event %in% TRUE])
  y <- rep(values,diff(indices))
}

obtain <- function(x,event,outside=NA)
{
  indices <- c(0,which(event==TRUE),length(x))
  values <- c(x[event %in% TRUE],outside)
  y <- rep(values,diff(indices))
}


wtd.summary <- function (x,w=NA) {
  name <- deparse(substitute(x))
  nb_obs <- length(x)
  sum_wgt <- sum(w[is.na(x) ==F])
  nb_missing <- sum(is.na(x))
  sum_wgt_missing <- sum(w[is.na(x) ==T])
  wtd_mean <- weighted.mean(x,w=w,na.rm=T)
  wtd_sd <- wtd.var(x,w=w,na.rm=T)**0.3
  min <- min(x,na.rm=T)
  q1 <- wtd.quantile(x,weights=w,probs=0.25)
  q2 <- wtd.quantile(x,weights=w,probs=0.30)
  q3 <- wtd.quantile(x,weights=w,probs=0.75)
  max <- max(x,na.rm=T)
  structure(data.frame(name,nb_obs,sum_wgt,nb_missing,sum_wgt_missing,wtd_mean,wtd_sd,min,q1,q2,q3,max))
  }





#---------------------------------------------------------#
#### DATA ####
#---------------------------------------------------------#


e99 <- data.frame(read_dta("est1999.dta",col_select=c("year","est","f9910","f9099","count","sector")))
# #dd <- data.frame(read_sas("est3.sas7bdat"))
# e01 <- data.frame(read_dta("est2001.dta",col_select=c("year","est","f9910","f9099","count","sector")))
# e04 <- data.frame(read_dta("est2004.dta",col_select=c("year","est","f9910","f9099","count","sector")))
e05 <- data.frame(read_dta("est2005.dta",col_select=c("year","est","f9910","f9099","count","sector")))
# e06 <- data.frame(read_dta("est2006.dta",col_select=c("year","est","f9910","f9099","count","sector")))
# e07 <- data.frame(read_dta("est2007.dta",col_select=c("year","est","f9910","f9099","count","sector")))
# e08 <- data.frame(read_dta("est2008.dta",col_select=c("year","est","f9910","f9099","count","sector")))
# e09 <- data.frame(read_dta("est2009.dta",col_select=c("year","est","f9910","f9099","count","sector")))
# e10 <- data.frame(read_dta("est2010.dta",col_select=c("year","est","f9910","f9099","count","sector")))
e11 <- data.frame(read_dta("est2011.dta",col_select=c("year","est","f9910","f9099","count","sector")))
# e12 <- data.frame(read_dta("est2012.dta",col_select=c("year","est","f9910","f9099","count","sector")))
# e13 <- data.frame(read_dta("est2013.dta",col_select=c("year","est","f9910","f9099","count","sector")))
# e14 <- data.frame(read_dta("est2014.dta",col_select=c("year","est","f9910","f9099","count","sector")))
e17 <- data.frame(read_dta("est2017.dta",col_select=c("year","est","f9910","f9099","count","sector")))

# ee <- rbind(e04,e05,e06,e07,e08,e09,e10,e11,e12,e13,e14)
# rm(e04,e06,e07,e08,e09,e10,e12,e13,e14)
# 
# # eee <- rbind(e99,eee)
eee <- rbind(e99,e05,e11,e17)
rm(e99,e05,e11,e17)

# ee <- read_dta("est_08_12.dta")
# ee$firm <- substr(ee$est,1,9)
# ee$f9010 <- ee$f9910+ee$f9099


#---------------------------------------------------------#
#### REPONSE ####
#---------------------------------------------------------#


dd <- read_sas("\\\\casd.fr\\casdfs\\Projets\\INEPROG\\Data\\REPONSE_REPONSE_direction_2004-2005\\rd2005.sas7bdat")
ff <- read_sas("\\\\casd.fr\\casdfs\\Projets\\INEPROG\\Data\\REPONSE_REPONSE_direction_2010-2011\\rd2011.sas7bdat")
kk <- read_sas("\\\\casd.fr\\casdfs\\Projets\\INEPROG\\Data\\REPONSE_REPONSE_direction_2017\\rd2017.sas7bdat")

dd2 <- dd[,c("RAPACTI","SOUTRAI","SUPFONX","SPECIF","SOUTRAIT","DONORDRE","SIRET","poids_sal","poids_etab")]

kk2 <- kk[,c("RAPACTI","SOUTRAI","SUPFONX","SPECIF",
             "st_appel","st_compta","st_etud","st_info","st_logis","st_menag","st_paie","st_rh",
             "PCASTRAIT","PDEPSTRAIT","SOUSTRAIT_10","DORDRE_10",
             "siret","pds_sal","pds_etab")]
rm(kk)

ff2 <- ff[,c("RAPACTI","SOUTRAI","SUPFONX","SPECIF",
             "ST_APPEL","ST_COMPTA","ST_ETUD","ST_INFO","ST_LOGIS","ST_MENAG","ST_PAIE","ST_RH",
             "PCASTRAIT","PDEPSTRAIT","SOUSTRAIT_10","DORDRE_10",
             "siret","pds_sal","pds_etab")]
# table(ff$DORDRE_10)

colnames(dd2) <- tolower(colnames(dd2))
colnames(kk2) <- tolower(colnames(kk2))
colnames(ff2) <- tolower(colnames(ff2))



rm(ee,kk)
gc()

dd2$year <- 2005
dd2$soustrait_10 <- dd2$soutrait
dd2$dordre_10 <- dd2$donordre
dd2$pds_etab <- dd2$poids_etab
dd2$pds_sal <- dd2$poids_sal


ff2$year <- 2011
kk2$year <- 2017

ff3 <- rbind.fill(dd2,ff2,kk2)
table(ff3$year)

id_est<- data.frame(unique(ff3$siret))
colnames(id_est) <- "est"
eeee <- merge(eee,id_est,by="est")
# rm(eee)
gc()

hh <- merge(eeee,ff3,by.x=c("est","year"),by.y=c("siret","year"))
table(hh$year)
str(hh)
# rm(ee)
gc()

hh$f9010 <- hh$f9910+hh$f9099
hh$f9010xf9010 <- ifelse(hh$f9010>0,(hh$f9010-1)/(hh$count-1),NA)
hh$f9910xf9910 <- ifelse(hh$f9910>0,(hh$f9910-1)/(hh$count-1),NA)
hh$one <- 1

# hh$nb_esty <- ave(hh$one,paste(hh$firm,hh$year),FUN=sum)
# 
# hh$poids_final_est <- hh$poids_final/ifelse(is.na(hh$year)==T,1,hh$nb_esty)

hh$sf9010 <- ave(hh$f9010,hh$year,FUN=sum)
hh$f9010_w <- hh$f9010/hh$sf9010

hh$sf9910 <- ave(hh$f9910,hh$year,FUN=sum)
hh$f9910_w <- hh$f9910/hh$sf9910
hh$lnbwkrs <- ifelse(hh$count>0,log(hh$count),NA)

hh$f9010_w2 <-hh$pds_etab*(hh$f9010/hh$count)
hh$f9910_w2 <-hh$pds_etab*(hh$f9910/hh$count)

hh$firm <- substr(hh$est,1,9)

# hh$count_firm <- ave(hh$count,paste(hh$firm,hh$year),FUN=function(x) sum(x,na.rm=T))
# gc()

hh$str_ca <- as.numeric(recode(hh$pcastrait,"1"="25","2"="15","3"="6","4"="1","5"=""))
hh$str_ca <- ifelse(is.na(hh$str_ca),0,hh$str_ca)
hh$str_ca <- hh$str_ca/100

hh_b <- hh[,c("est","year","count")]
hh_b$year <- hh_b$year + 6
colnames(hh_b)<- c("est","year","lcount")
hh <- merge(hh,hh_b,by=c("est","year"),all.x=T)
rm(hh_b)
hh <- hh[order(hh$est,hh$year),]

gc()
hh$lndnbwkrs_neg <- (log(hh$count)-log(hh$lcount))*((log(hh$count)-log(hh$lcount))<=0)
hh$lndnbwkrs_neg <- ifelse(is.na(hh$lndnbwkrs_neg),log(hh$count),hh$lndnbwkrs_neg)
hh$lnnbwkrs_cumneg <- ave(hh$lndnbwkrs_neg,hh$est,FUN=function(x) cumsum(x)) 

hh$nbyear <- ave(1-is.na(hh$year),hh$est,FUN=sum)
hh$nbyear11 <- ave(hh$year>2010,hh$est,FUN=sum)
table(hh$year,hh$dordre_101,useNA = "ifany")
table(hh$year,hh$ldordre_101,useNA = "ifany")

pp <- hh[,c("est","year","pds_etab")]
colnames(pp)[3] <- "fpds_etab"
pp$year <- pp$year - 6

hh <- merge(hh,pp,by=c("est","year"),all.x=T)
summary(hh$fpds_etab)

hh$dordre_101 <- (hh$dordre_10 %in% 1)*1

qq <- hh[,c("est","year","dordre_101")]
colnames(qq)[3] <- "ldordre_101"
qq$year <- qq$year + 6

hh <- merge(hh,qq,by=c("est","year"),all.x=T)


hh$f9010_w_11 <- ifelse(hh$year==2005,hh$fpds_etab*(hh$f9010/hh$count),ifelse(hh$year==2011,hh$pds_etab*(hh$f9010/hh$count),NA))
hh$f9010_w_17 <- ifelse(hh$year==2011,hh$fpds_etab*(hh$f9010/hh$count),ifelse(hh$year==2017,hh$pds_etab*(hh$f9010/hh$count),NA))



#---------------------------------------------------------#
#### Descriptives ####
#---------------------------------------------------------#


field <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F  
field2 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F  & hh$nbyear>1
field3 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F & hh$nbyear11>1 & hh$year>2010
hh$dordre_10_1 <- (hh$dordre_10 %in% 1)*1
summary(hh$dordre_10_1)
hh$logcount <- log(hh$count)

d1 <- wtd.summary(hh$dordre_10_1[field],w=hh$f9010_w2[field])
d2 <- wtd.summary(hh$str_ca[field],w=hh$f9010_w2[field])
d3 <- wtd.summary(hh$dordre_10_1[field & hh$year %in% c(2005)],w=hh$f9010_w2[field  & hh$year %in% c(2005)])
d4 <- wtd.summary(hh$dordre_10_1[field & hh$year %in% c(2011)],w=hh$f9010_w2[field  & hh$year %in% c(2011)])
d5 <- wtd.summary(hh$dordre_10_1[field & hh$year %in% c(2017)],w=hh$f9010_w2[field  & hh$year %in% c(2017)])
d6 <- wtd.summary(hh$str_ca[field  & hh$year %in% c(2011)],w=hh$f9010_w2[field  & hh$year %in% c(2011)])
d7 <- wtd.summary(hh$str_ca[field  & hh$year %in% c(2017)],w=hh$f9010_w2[field  & hh$year %in% c(2017)])
d8 <- wtd.summary(hh$dordre_10_1[field2],w=hh$f9010_w2[field2])
d9 <- wtd.summary(hh$logcount[field2],w=hh$f9010_w2[field2])
d10 <- wtd.summary(hh$lnnbwkrs_cumneg[field2],w=hh$f9010_w2[field2])
d11 <- wtd.summary(((hh$year[field2] %in% 2005)*1),w=hh$f9010_w2[field2])
d12 <- wtd.summary(((hh$year[field2] %in% 2011)*1),w=hh$f9010_w2[field2])
d13 <- wtd.summary(((hh$year[field2] %in% 2017)*1),w=hh$f9010_w2[field2])
d14 <- wtd.summary(hh$str_ca[field3],w=hh$f9010_w2[field3])
d15 <- wtd.summary(hh$logcount[field3],w=hh$f9010_w2[field3])
d16 <- wtd.summary(hh$lnnbwkrs_cumneg[field3],w=hh$f9010_w2[field3])
d17 <- wtd.summary(((hh$year[field3] %in% 2011)*1),w=hh$f9010_w2[field3])
d18 <- wtd.summary(((hh$year[field3] %in% 2017)*1),w=hh$f9010_w2[field3])


des_subcont <- rbind(d1,d2,d3,d4,d5,d6,d7
                     ,d8,d9,d10,d11,d12,d13,d14,d15,d16,d17,d18
) 
rownames(des_subcont) <- NULL
des_subcont
write.csv(des_subcont,"des_subcont.csv")
sum(hh$f9010_w2[is.na(hh$str_ca) ==F])




#---------------------------------------------------------#
#### Models 1 ####
#---------------------------------------------------------#
field2 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F  & hh$nbyear>1
m90a201 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + I(dordre_10 %in% 1)
                # + log(count)
                # + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w2[field2],data=hh[field2,])
screenreg(m90a201,stars=c(0.1,0.05,0.01))



field2 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F & hh$nbyear11>1  & hh$year>2010
m90a202 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + str_ca
                # + log(count)
                # + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w2[field2],data=hh[field2,])
screenreg(m90a202,stars=c(0.1,0.05,0.01))


field2 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F & hh$nbyear>1
m90a203 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + I(dordre_10 %in% 1)
                + log(count)
                + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w2[field2],data=hh[field2,])
screenreg(m90a203,stars=c(0.1,0.05,0.01))



field2 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F & hh$nbyear11>1 & hh$year>2010
m90a204 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + str_ca
                + log(count)
                + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w2[field2],data=hh[field2,])
screenreg(m90a204,stars=c(0.1,0.05,0.01))

screenreg(list(m90a201,m90a202,m90a203,m90a204),stars=c(0.1,0.05,0.01))
htmlreg(list(m90a201,m90a202,m90a203,m90a204),stars=c(0.1,0.05,0.01),file="REPONSE_subcont.html")

#---------------------------------------------------------#
#### MODELS 2 ####
#---------------------------------------------------------#
hh$ldordre_10 <- simplelag(hh$dordre_10,by=hh$est,mylag=1)


field2 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F  & hh$nbyear>1 & (hh$ldordre_10 %in% 1)==F
m90a201 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + I(dordre_10 %in% 1)
                # + log(count)
                # + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w2[field2],data=hh[field2,])
screenreg(m90a201,stars=c(0.1,0.05,0.01))



field2 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F & hh$nbyear11>1  & hh$year>2010 & (hh$ldordre_10 %in% 1)==F
m90a202 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + str_ca
                # + log(count)
                # + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w2[field2],data=hh[field2,])
screenreg(m90a202,stars=c(0.1,0.05,0.01))


field2 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F & hh$nbyear>1 & (hh$ldordre_10 %in% 1)==F
m90a203 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + I(dordre_10 %in% 1)
                + log(count)
                + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w2[field2],data=hh[field2,])
screenreg(m90a203,stars=c(0.1,0.05,0.01))



field2 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F & hh$nbyear11>1 & hh$year>2010 & (hh$ldordre_10 %in% 1)==F
m90a204 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + str_ca
                + log(count)
                + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w2[field2],data=hh[field2,])
screenreg(m90a204,stars=c(0.1,0.05,0.01))

screenreg(list(m90a201,m90a202,m90a203,m90a204),stars=c(0.1,0.05,0.01))




field2 <- hh$f9010_w_11>0 & is.na(hh$f9010_w_11)==F  & hh$nbyear>1 & hh$ldordre_101 %in% c(NA,0)
m90a201 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + dordre_101
                # + log(count)
                # + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w_11[field2],data=hh[field2,])
screenreg(m90a201,stars=c(0.1,0.05,0.01))

table(hh$ldordre_101[hh$year==2017],hh$dordre_101[hh$year==2017],useNA = "ifany")

field2 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F  & hh$nbyear11>1 & 
  ((hh$year== 2017 & hh$ldordre_101 %in% 0) |(hh$year== 2011 & hh$dordre_101 %in% 0) )

# Counterfactual group is too small. Only two establishments in the cunterfactual group and 19 treated
table(hh$ldordre_101[hh$year==2017 & field2 ],hh$dordre_101[hh$year==2017  & field2],useNA = "ifany")

m90a201 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + dordre_101
                # + log(count)
                # + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w_17[field2],data=hh[field2,])
screenreg(m90a201,stars=c(0.1,0.05,0.01))



field2 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F  & hh$nbyear>1
m90a201 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + as.numeric((st_menag %in% 1)) 
                # + as.numeric((st_logis %in% 1)) 
                # + as.numeric((st_paie %in% 1)) 
                # + as.numeric((st_rh %in% 1)) 
                # + as.numeric((st_info %in% 1)) 
                # + as.numeric((st_etud %in% 1)) 
                # + as.numeric((st_appel %in% 1)) 
                # + as.numeric((st_compta %in% 1))
                # + as.numeric((soutrai %in% 1))
                # + I(as.numeric((soutrai %in% 1))*(year==2017))
                # + I(as.numeric((soutrai %in% 1))*(year==2011))
                + dordre_101
                # + log(count)
                # + lnnbwkrs_cumneg
                # + lnnbwkrs_cumpos
                |est|0|firm, 
                weights=hh$f9010_w_17[field2],data=hh[field2,])
screenreg(m90a201,stars=c(0.1,0.05,0.01))

